Abundance, biomass and species richness of macrozoobenthos along an intertidal elevation gradient

Abstract Ecology aims to comprehend species distribution and its interaction with environmental factors, from global to local scales. While global environmental changes affect marine biodiversity, understanding the drivers at smaller scales remains crucial. Tidal flats can be found on most of the world's coastlines and are particularly vulnerable to anthropogenic disturbances. They are important transient ecosystems between terrestrial and marine ecosystems, and their biodiversity provides important ecosystem services. Owing to this unique, terrestrial–marine transition, strong environmental gradients of elevation, sediment composition and food availability prevail. Here, we investigated which regional and local environmental factors drive the spatio‐temporal dynamics of macrozoobenthos communities on back‐barrier tidal flats in the East Frisian Wadden Sea. On the regional level, we found that species composition changed significantly from west to east on the East Frisian islands and that total abundance and species richness decreased from west to east. On the local abiotic level, we found that macrozoobenthos biomass decreased with higher elevation towards the salt marsh and that the total abundance of organisms in the sediment significantly increased with increasing mud content, while biodiversity and biomass were not changing significantly. In contrast to expectations, increasing Chl a availability as a measure of primary productivity did not result in changes in abundance, biomass or biodiversity, but extremely high total organic carbon (TOC) content was associated with a decrease in biomass and biodiversity. In conclusion, we found regional and local relationships that are similar to those observed in previous studies on macrozoobenthos in the Wadden Sea. Macrozoobenthos biomass, abundance and biodiversity are interrelated in a complex way with the physical, abiotic and biotic processes in and above the sediment.


| INTRODUC TI ON
One of the main aims in ecology is to understand the distribution of species and how they are influenced by abiotic and biotic environmental factors (Sutherland et al., 2013).Global environmental change has far-reaching consequences for marine biodiversity and ecosystem functioning (Bulling et al., 2010;Dulvy et al., 2003;Worm et al., 2006), but predicting diversity and community composition at a global scale might not help to understand factors that are driving the diversity and community composition at smaller scales (Willis & Whittaker, 2002).Responses of marine macrozoobenthos communities to multiple drivers are far from being understood (Cheung et al., 2009;Hoegh-Guldberg & Bruno, 2010;Ysebaert & Herman, 2002).At the regional scale, dispersal dynamics can play an important role in shaping the community composition and community properties (Palmer et al., 1996).On the local scale, environmental filtering, the prevention of establishment or persistence of an organism at a certain location can shape the community properties (Kraft et al., 2015).Furthermore, local food resource availability can also act as a filter by leading to the death of species (i.e., those unable to tolerate low food availability; Bauer et al., 2022).
Additionally, biotic interactions such as competition and predation and positive trophic and non-trophic interactions affect the communities (Bertness & Leonard, 2014;Ysebaert & Herman, 2002).Biotic interactions are changing with increasing temperatures, anthropogenic nutrient inputs, range expansion of or appearance of invasive species and exploitation by fishery (Beukema & Dekker, 2020).
Deepening our understanding of the interplay of drivers affecting species richness, community composition, species traits and ecological functions is crucial to mitigate global change consequences on marine ecosystems by identifying structuring environmental factors across various spatial and temporal scales.
Coastal areas, including the Wadden Sea, are dynamic and unique ecosystems bridging the terrestrial and open sea ecosystems.
Besides providing valuable economic (e.g., fish) and recreational ecosystem services (tourism), the large biomass of associated invertebrates supports the coastal food chain (de la Vega et al., 2018;Heip et al., 1995;Murray et al., 2019;Schückel et al., 2015;Sijtsma et al., 2019).The Wadden Sea comprises thousands of square meters of sand and mud flats located between barrier islands and the mainland coast (Essink et al., 1998;Flemming, 2012).Predatory fish and hundreds of thousands of migratory birds depend on macrozoobenthos as a food source and utilize this habitat as a nursery ground each year (Horn et al., 2020).Since 2009, the Wadden Sea has been declared a World Heritage Site due to its unique nature.Increasing habitat alteration (land reclamation and diking), exploitation (bottom trawling and shellfish dredging) and pollution (nutrients and chemicals) have changed the abundance and distribution of species in the Wadden Sea decreasing ecosystem functioning and the complexity of food webs (Eriksson et al., 2010;Lotze, 2005;Lotze et al., 2006).
Furthermore, extreme weather events, expected to increase in frequency in the future, such as warm or unusually cold winters (ice-winters), can lead to reduced recruitment and changes in the abundance of benthic species like bivalves and polychaetes (Bartels-Hardege & Zeeck, 1990;Beukema & Dekker, 2005).
Mostly east-ward directed winds in the German Bight (Figure 1) induce the residual circulation to be cyclonic, going from west to east in the southern parts of the Wadden Sea (Otto et al., 1990;Staneva et al., 2009).During flood, the circulation is directed eastwards, changing to a westward circulation during ebb where water from the tidal flats is sucked out via the tidal channels (Meyerjürgens et al., 2019;Stanev et al., 2007).Hydrodynamic models and riverine nutrient inputs suggest a decreasing gradient of nutrient levels and possibly also larval supply from west to east along the Frisian islands, from Norderney over Spiekeroog to Wangerooge (Respondek et al., 2022;Stanev et al., 2003;van Beusekom et al., 2009van Beusekom et al., , 2019)), which could have an influence on the macrozoobenthos community composition at the islands.
Benthic macrozoobenthos is well adapted to the dynamics and harsh conditions of the tidal flats that are characterized by diurnal tides, and therefore constantly changing water levels, which leads to strong environmental gradients (Le Hir et al., 2000;Schückel et al., 2013;Sprung et al., 2001).Tidal flat organisms exhibit distinct distribution patterns within specific sections of the tidal flat, influenced by environmental gradients, primarily elevation, sediment composition and total organic carbon (Beukema & Dekker, 2009;Compton et al., 2013;Reise, 2002;Schückel et al., 2013;Ysebaert et al., 2003).Strong abiotic gradients induce a typical zonation of tidal flats with low species richness and high abundance of small organisms at the high-tide line, close to the salt marsh (Beukema & Dekker, 2009).With decreasing elevation, species richness as well as biomass increase and then decline again towards exposed areas in the subtidal (Beukema & Dekker, 2009;Donadi et al., 2015;Schückel et al., 2013).Tidal flats can be classified into mud-, mixedand sand flats based on the content of mud or sand in the sediment (Flemming, 2000;Reineck & Siefert, 1980).Sediment composition can play a major role in structuring macrozoobenthos communities since many species are associated with a specific type of sediment, for example, for building tubes but also by reflecting the hydrodynamic regime and food availability (Armonies, 2021;Kröncke, 2006;Nehmer & Kröncke, 2003;Reise, 2002;Reiss & Kröncke, 2001;Schückel et al., 2013;Schückel & Kröncke, 2013).Overall, macrozoobenthos communities are predominantly influenced by two environmental factors: elevation and sediment.
Biotic conditions, such as primary production (Chl a biomass), can have influence on the macrobenthic communities by being a

T A X O N O M Y C L A S S I F I C A T I O N
Biodiversity ecology, Community ecology, Ecosystem ecology, Soil ecology, Spatial ecology, Taxonomy, Trophic interactions primary food source (Beukema & Cadée, 1997;Puls et al., 2012;Schückel et al., 2013).Next to seasonal variation, human activities, like increased nutrient run-off from rivers, changes in hydrodynamics due to diking, loss of seagrass and bivalve fisheries, also influence primary production (Eriksson et al., 2010;Philippart et al., 2007).The Wadden Sea undergoes eutrophication and oligotrophication cycles (Burson et al., 2016;Lenhart et al., 2010;van Beusekom et al., 2019;van Raaphorst & de Jonge, 2004), affecting the total biomass and community composition of primary producers, which, in turn, cascade to higher trophic levels (Beukema & Dekker, 2022;Cloern, 2001;van Roomen et al., 2012).Nutrient availability is a central factor for the occurrence of microphytobenthos, consisting of unicellular eukaryotic algae and cyanobacteria (Hope et al., 2020;MacIntyre et al., 1996), and consequently, for the functioning of coastal soft-sediment ecosystems.
Microphytobenthos contributes to a large extent to the primary production of tidal flats and supports higher trophic levels, like macrozoobenthos, with essential fatty acids and energy by contributing large fractions to their nutrition (Christianen et al., 2017;Hope et al., 2020;MacIntyre et al., 1996).Microphytobenthos biomass generally increases with elevation and cohesiveness of the sediment, most likely due to the longer exposure time and therefore light availability, higher nutrient concentration and less resuspension of the sediment (Colijn & Dijkema, 1981;Underwood & Kromkamp, 1999;Yong et al., 2022).While opportunistic species do best in habitats with high organic enrichment, equilibrium species can better cope with lower concentrations (Linton & Taghon, 2000;Pearson & Rosenberg, 1978).Species performance in habitats with different levels of enrichment can depend on their physiological and behavioural adaptations (Linton & Taghon, 2000;Pearson & Rosenberg, 1978).
Despite the importance of understanding the distribution of macrozoobenthos on different spatial scales and under varying environmental conditions, data availability of the community composition of macrozoobenthos in combination with environmental and functional parameters (e.g., biomass) from the East Frisian islands is limited (Bergfeld, 1999;Dörjes et al., 1986;Hodapp et al., 2014;Kröncke, 1996;Nehmer & Kröncke, 2003;Reiss & Kröncke, 2001;Singer et al., 2023), and an evaluation of major environmental drivers across different islands is mainly missing.
The objective of this study is to quantify the abundance, biomass and diversity of benthic macrozoobenthos at both regional and local spatial levels, as well as in two seasons, in relation to abiotic and biotic environmental factors on tidal flats in the German Wadden Sea.Thus, we formulated the hypotheses as follows: • H1: At the regional scale, the influence of dispersal limitation leads to fewer shared species between eastern and western communities, resulting in pronounced dissimilarity in their community composition.
• H2.Also, at the regional scale, nutrient-induced species sorting, leads to a gradient in macrozoobenthos abundance, biomass and biodiversity, with higher levels observed in western areas compared with eastern areas.
F I G U R E 1 Map shows the three East Frisian islands, Norderney, Spiekeroog and Wangerooge with three transects (A = yellow, B = red and C = blue) at each island and three sample stations within each transect.At each sampling station, three samples were taken in spring and summer.The map was created with QGIS (version 3.24.1)and map data from OpenStreetMap.
• H4: At the local scale, biotic filtering in terms of resource availability, displayed by microphytobenthos (Chl a) and TOC, causes an increase in abundance, biomass and diversity with increasing resource availability (Pearson & Rosenberg, 1978;Schückel et al., 2013).However, at extreme food availability, species richness and diversity decline, followed by a decrease in biomass due to an increase in abundant but small-sized species (Pearson & Rosenberg, 1978).Additionally, recruitment of larvae in spring occurs at sites that provide suitable abiotic and biotic conditions, resulting in higher abundance and biomass in summer (Beukema, 1974).We have no particular expectations relative to diversity and season.
The islands are part of the Wadden Sea National Park of Lower Saxony, located in the Southern North Sea of Germany.The tidal flats, which are on the southern side of the barrier islands, facing the mainland, are shaped by semidiurnal tides and relatively low current velocities (Hild et al., 1999).

| Macrozoobenthos community
To cover a broad range of environmental conditions (Table A2), we sampled a total of nine transects with 27 sampling stations, that is, three transects perpendicular to island coast at each island comprised of three sampling stations each.Each sampling station was marked with a DGPS to estimate the elevation of the study plot.
The sampling stations covered a depth range from the high-tide line to the middle-tide line (Figure 1).To cover the spatial and temporal changes in the macrozoobenthos community, sampling took place once in spring and once in summer of 2019 (Table A1).At each of the sampling stations, we took three subsamples for macrozoobenthos with hand corers that had a diameter of 10 cm (0.0078 m 2 ) and a height of 20 cm, which were driven into the sediment at low tide.
The three cores were placed in separate plastic bags, transported to the shore and washed over a 500 μm sieve.Remaining organisms were then placed into separate plastic containers and fixated using 4% buffered-formaldehyde solution for later identification.In the laboratory, samples were stained with rose bengal and identified to species level (or the lowest possible taxonomic unit) under a stereomicroscope (Leica MZ 95; Leitz, Wetzlar, Germany).For smaller organisms, an optical microscope (Leica DM LS 2) was used.
After taxonomic identification, abundances of each species and wet biomass were determined.The World Register of Marine Species (WoRMS) taxon match tool was used to standardize the nomenclature (http:// www.marin espec ies.org/ aphia.php?p= match ).Not all organisms could be determined to species level; this was especially the case for oligochaetes, insect larvae and juvenile specimens.
Therefore, species richness calculations were partly based on number of higher taxonomic units, resulting in a conservative estimation of diversity.Although only 15.2% of the specimens were identified to species level, we decided to still use the term species richness in data evaluation and discussion.

| Total organic carbon and sediment composition
To measure the total organic carbon (TOC) content, one surface sediment sample (<2 cm depth) per station was taken and frozen at −20°C until further analysis.The samples were then freeze-dried in a Lyovac GT2 (STERIS GmbH, Huerth, Germany) and subsequently ground using a planetary ball mill (PM200; Retsch GmbH, Haan, Germany).Twenty-five to seventy-five milligrams of the material was weighed in silver cups and placed in a steaming 37%-HCl desiccator for 24 h to remove carbonates.After that, one drop of 10%-HCl was added to remove all carbonate leftovers.For evaporating the leftover HCl, the samples were placed in an oven at 50°C for 12 h.Subsequently, the silver cups were folded into small packages and analysed for total organic carbon in a CHNS-Elemental-Analyzer 'vario EL cube' using Loess as a standard (ELEMENTAR Analysensysteme GmbH Heraeus, Langenselbold, Germany).
To analyse the particle size distribution of the sediment at each station, a sediment sample of the upper 5 cm was taken with a hand corer (diameter 10 cm), placed in a plastic bag and frozen at −20°C until further analysis.Samples were thawed and mixed with a spatula and 20-35 g of sediment were dissolved in a 4%-sodiumpolyphosphate solution and sieved through a 2 mm sieve to remove large particles.To disperse the clay, the mixture was shaken overnight and subsequently measured with a laser diffraction spectrometer (HORIBA LA-950).Sediment fraction of mud (<0.063 mm) was obtained in volume per cent.

| Inorganic nutrients and primary production
We used the dissolved inorganic nutrients and Chl a biomass informed by Yong et al. (2022), who took samples at the exact same locations and time as our macrozoobenthos samples.Dissolved inorganic nutrients were sampled using rhizon samplers (CSS 5 cm; 0.15 μm membrane pore size; 19.21.23F;Rhizosphere Research Products, Wageningen, Netherlands) that were inserted 2 cm in the sediment and connected to syringes.Samples were analysed for reactive dissolved phosphate (DIP) and dissolved inorganic nitrogen (ammonium, nitrite and nitrate combined) using a Scalar analytical auto-analyzer (San++ System, Scalar Analytical, Breda, The Netherlands).
For Chl a determination, three sediment cores were taken at each sampling station from the upper 2 cm of the sediment with a cut-off syringe (diameter 26 mm) and pooled together (total volume 31.85 cm 3 ).In the field, samples were stored in the dark, cooled with ice packs and immediately frozen at −80°C after returning to the laboratory.The extraction and analysis of the pigment Chl a was performed after the method of Thrane et al. (2015) with slight modifications that can be found in Yong et al. (2022).In summary, Chl a was extracted from sediment subsamples (5 g) in 96% ethanol and the resulting supernatants were measured in triplicates using a Synergy MX plate reader (BioTek Instruments, Vermont, USA).Spectral scans were carried out between 300 and 800 nm with a resolution of 1 nm to measure Chl a concentrations, which is a proxy for total algal biomass.Results from the plate reader were analysed in R (R Core Team, 2021) with the R script provided by Thrane et al. (2015) that is based on a modified Gauss-peak spectra (GPS) method.

| Data analysis
Univariate and multivariate data were analysed using the program R (R Core Team, 2021;version 4.0.3 and RStudio Team, 2022).
Abundance data of macrozoobenthos in each subsample were standardized by total abundance of the subsample, because multivariate analysis can be sensitive to absolute abundances, and a Bray-Curtis dissimilarity matrix was calculated.To test the differences in community composition between islands (H1), we performed a PERMANOVA using the adonis2 function from the 'vegan' package with 9999 permutations (Oksanen et al., 2007).
To visualize the similarity or dissimilarity of the macrozoobenthos communities, a non-metric multidimensional scaling (NMDS) was performed.The resulting two-dimensional ordination plot displays the samples sorted relative to their dissimilarity, with dissimilar samples further apart and similar samples closer together.For the NMDS, the function metaMDS from the 'vegan' package with 9999 permutations was used (Oksanen et al., 2007).In order to fit the environmental data onto the ordination, we used the envfit function from the 'vegan' package with 9999 permutations and rm.na set to true, to calculate the correlation of environmental variables with the ordination axes.Visualization was realized with ggplot2 (Wickham et al., 2016), and all significant environmental variables are shown as vectors in Figure 2, Table 1.We calculated species richness and the inverse Simpson diversity index using the function hill_taxa from the 'hillR' package (Chao et al., 2014;    , 2018).Species richness (Hill number, q = 0), the inverse of Simpson's diversity index (Hill number, q = 2), the total number of individuals and the biomass were calculated per subsample.Data exploration was carried out following the protocol described by Zuur et al. (2010).Species richness and the inverse Simpson index were transformed logarithmically before analysis to ensure normality, while biomass was ln (x + 1) transformed.All fixed effects were scaled prior to analysis because fixed effects were on different scales.
We fitted separate models for the response variables abundance, biomass, species richness and diversity using linear mixed effect models (LMMs) for normal distributed response variables (diversity and species richness) and generalized mixed effect models (GLMMs) for non-normal distributed response variables (abundance and biomass) to test the hypothesis H2, H3 and H4 in combination using longitude (continuous), latitude (continuous), elevation (continuous), mud content (volume per cent), season (categorical), TOC (continuous) and Chl a (continuous) as explanatory variables.As link functions, we chose nbinom2, which is a negative binomial distribution with a quadratic parametrization, for the abundance (Hardin & Hilbe, 2018) and gaussian for species richness and inverse Simpson index.For biomass, which includes continuous but non-negative values, we chose a tweedie distribution that was shown to be often appropriate for this kind of data (Foster & Bravington, 2013;Niku et al., 2017).For LMMs, the 'lme4' package (Bates et al., 2015) was used; for GLMMs, the 'glmmTMB' package was used (Brooks et al., 2017).Explanatory variables were checked for multicollinearity using variance inflation factors (VIF) obtained by the vif function from the 'car' package (Fox et al., 2012).VIFs indicated a high correlation of latitude and longitude.Therefore, we removed latitude as explanatory variable from the model.All remaining VIFs were lower than the threshold of 3 (Zuur et al., 2010).To make models converge, we used the 'bobyqa' optimizer which allows for a series of convergence checks.We increased the number of iterations to 200,000 using the 'maxfun' argument in the LMMs to help achieve convergence when the default value is not sufficient.In GLMMs, we used the optimizer 'nlminb' and increased iterations to 200,000.Marginal and conditional R 2 were calculated based on Johnson (2014), for assessing the goodness-of-fit.
To incorporate the dependency among observations from the same island and sampling station, we included island and sampling station as a random effect.However, the models for abundance, biomass and species richness as response variables resulted in singular fits, caused by the lack of spatial variablesassociated variance, and therefore, we simplified the random effect structure down to sampling station alone for these models (Barr et al., 2013).All models were checked for over-and underdispersion, zero-inflation, homogeneity of variance and also the suitability of chosen error distributions using the DHARMa package (Hartig, 2022).To display the most important information of the models, the package 'sjPlot' was used (Lüdecke, 2022).To visualize the outcomes of the lmer and glmmTMB models, the predictions of the models were extracted using the 'sjPlot' package and combined with the raw data using the 'ggplot2' package (Wickham et al., 2016).

| Spatial and temporal variability of macrozoobenthos
We collected and identified a total of 31,908 organisms comprising 45 different taxa overall (Table A3).Twenty-five species of polychaetes, seven crustaceans, seven oligochaetes, five bivalves, two gastropods and one insect species were found during the study period.The species richness in subsamples varied between 5 and 18 (mean ± SD: 11.2 ± 2.84, n = 162).

| Regional differences in the macrozoobenthos community composition (H1)
Macrozoobenthos communities differed significantly between islands (Figure 2), but only a small fraction of the variability was explained by the factor 'island' (PERMANOVA, p < .0001,R 2 = .129,Table 1).Norderney and Wangerooge communities clustered further apart from each other, while the community of Spiekeroog was in-between (Figure 2

| Effect of regional-scale species sorting on abundance, biomass and diversity (H2)
Total abundance (Figure 3a, p = .002)and species richness (Figure 3c, p < .001)decreased significantly from west to east, while biomass showed a slight significant increase (Figure 3b, p = .014)in the opposite direction.Diversity showed no significant change from west to east (p = .121).
Total abundance showed no change with elevation (p = .986,Table 2).
Sediment mud content had a significant positive effect on the total abundance (Figure 3j, p = .012,Table 2).Biomass, species richness and diversity depicted positive but non-significant changes with mud content in the sediment (p = .43,p = .273,p = .16,Table 2).

| Effects of local biotic factors and season on the abundance, biomass, species richness and diversity of macrozoobenthos (H4)
Algal biomass (Chl a) had a non-significant positive effect on total abundance and biomass (p = .351,p = .068,Table 2) and a non-significant negative effect on species richness and diversity (p = .195,p = .161,Table 2).The amount of TOC in the sediment had a significant negative effect on the biomass (Figure 3g, p = .008,Table 2), species richness and diversity (Figure 3h,i, p = .046,p = .007,Table 2).
Seasonal changes were observed in total abundance, which increased significantly towards summer (Figure 3k, p < .001,Table 2), while at the same time, biomass showed a positive but non-significant change towards summer (p = .068,Table 2).In contrast, the diversity declined significantly towards summer (Figure 3l,     takes both the fixed and random effects into account.Similar to the R 2 , the intra-class correlation coefficients (ICCs) give information on how much of the proportion of variance is explained by a grouping (random) factor in multilevel/hierarchical data (Nakagawa et al., 2017).The random effect variance σ 2 stands for the mean random effect variance of the model.The random intercept variance, or between-subject variance (τ 00 ), indicates how much groups differ from each other.N stands for the number of random effect groups.The family nbinom2 is a negative binomial distribution with quadratic parameterization (Hardin & Hilbe, 2018).***p < .001;**p < .01;*p < .05. p = .003,Table 2) and species richness did not change significantly (p = .969,Table 2) with season.
Marginal R 2 values indicate how much of the variance is explained by the fixed effects while the conditional R 2 values explain the variance of both the fixed and the random effects.The marginal R 2 values in the models ranged between .187 and .496and between .574 and .697for conditional R 2 values (Table 2), indicating that the random effects also explained some of the variation in the models.

| DISCUSS ION
We investigated how community composition, abundance, biomass, species richness and diversity changed on the regional level from west to east and how local environmental factors (elevation, mud content, Chl a and TOC) and season affected these community properties.We found that (H1) species composition differed significantly between islands and (H2) that abundance and species richness decreased from west to east (H3).Biomass decreased with elevation, while diversity and species richness increased with elevation, while abundance increased with mud content (H4).Chl a showed no significant effects, but TOC affected biomass, species richness and diversity negatively.Total abundance increased from spring to summer, but inverse Simpson diversity decreased.

| Regional-scale dispersal limitation and nutrient-induced species sorting: Implications for macrozoobenthos community composition and characteristics (H1)
Macrozoobenthos species composition might be shaped by dispersal limitations between the islands, resulting in significantly different species compositions at the islands and with western communities being less similar to more eastern communities (Figure 2, Table 1).Species composition of macrozoobenthos can be influenced by larval dispersal and larval settlement (Chust et al., 2016).Most macrozoobenthos species disperse via planktonic larvae but can also drift from one location to the other as juveniles or adults (Armonies, 1999).Therefore, the anticlockwise residual currents, the tides, the position of the tidal channels and the topography in the region (Stanev et al., 2003) could be one factor influencing the dispersal and recruitment of larvae (Bouma et al., 2001) and hence the community composition on the tidal flats of the islands.Another factor, known to affect macrozoobenthos community composition, is anthropogenic influences at the islands.After severe storm surges in the past, the East Frisian is- Sea that were explained mainly by mean grain size, microphytobenthos biomass and exposure time.Overall, functional traits, such as recruitment success, tolerance to abiotic conditions like the sediment characteristics (Anderson, 2008;Compton et al., 2009) or inundation times (Kraan et al., 2010) but also the local species pool and the food availability, can change the community composition (Beukema & Dekker, 2005;Schückel et al., 2013).

| Regional changes of abundance, species richness and biomass (H2)
We also found decreasing abundance and species richness from west to east.Abundance and species richness is like the community composition, influenced by dispersal (Mouquet & Loreau, 2003).
Due to the prevailing hydrography, it is likely that the recruitment of different species is higher in the west compared with the east, thus increasing the abundance and diversity (Stanev et al., 2003).Even though the inverse Simpson index was not significantly affected by longitude, the direction of the trend was also negative (Table 2).
Other than expected in hypothesis H2, biomass increased slightly towards the east.Compton et al. (2013) also observed increasing biomass from west to east in the Dutch Wadden Sea and attributed it to increasing microphytobenthos in the east, which is not the case in our study (Figure A1F).One possible explanation is the increase in suitable habitats for species dominating the biomass on tidal flats, such as the common cockle, C. edule, that prefers sandier sites (Figure 2; Kraan et al., 2010).

| Impact of local abiotic drivers (mud content and elevation) on macrozoobenthos: Abundance, biomass and biodiversity patterns (H3)
The decrease in biomass with elevation partly supported H3.
One of the reasons, for the biomass decline might be due to the higher abundance of larger bivalves at mean-tide level (Beukema & Dekker, 2009), often contributing in large proportions to the biomass of tidal flats (Daggers et al., 2020;Essink et al., 1998).Often, filter feeders like C. edule profit from longer inundation times at lower elevation levels having more time to filter phytoplankton, as well as surface (deposit) feeders, which are mostly restricted to times of inundation to feed (Beukema & Dekker, 2009;Van Colen et al., 2010).
A relatively small number of species were responsible for a similar pattern in a study from Hodapp et al. (2014)  We observed significant positive changes in species richness and diversity with increasing elevation.This pattern contrasted with the traditional reported zonation of tidal flats where diversity is low near the salt marsh and increases towards low tide levels (Beukema & Dekker, 2009).Species at the border of the salt marsh tend to be characterized by small body size.Many studies investigating the zonation of macrozoobenthos, used a mesh size of 1 mm (e.g., Beukema & Dekker, 2009), while a mesh size of 0.5 mm was used in the current study.The study of Bachelet (1990) compared the efficacy of different mesh sizes and found that only 15% of polychaetes were captured in studies using 1 mm mesh size.Sieves with 0.5 mm mesh, like in this study, resulted in 53% of capture of polychaetes (Bachelet, 1990;Reise et al., 1994).This could result in different patterns of species richness and diversity in intertidal flats.
The transition between tidal flat and salt marsh, is also a transition zone between two ecological communities called ecotone (Bearup & Blasius, 2017).Ecotones are often areas were species of two different communities overlap, which can lead to higher biodiversity compared with the adjacent habitats (Kark & van Rensburg, 2006), which could also be the case here.
The tidal flats on Wangerooge were characterized as sandy (mud content <10%), while the other two islands have mixed sediment types (mud content 10-50%; Flemming & Davis Jr., 1994).Even though mud content did not explain total biomass, species richness and diversity as we expected in H4, it had a significant positive effect on the total abundance in our study (Figure 3j, Table 2).Therefore, H4 can be partly accepted, as abundance increased with increasing mud content.This indicates a shift to a higher abundance of smaller species, which was also demonstrated by the dominance of oligochaetes and small polychaetes like Baltidrilus costatus and Manayunkia aestuarina at sites with high mud content.This pattern is particularly strong for Wangerooge, where a high sand content and the lack of the typical zonation with higher mud content at the shore caused by higher hydrodynamic stress resulted in lower macrozoobenthos abundance and species richness (Friedrichs, 2012;Schückel et al., 2013;van der Wal, Ysebaert, & Herman, 2017).Overall, human modifications of the coast caused current velocities to increase, which has led to a depletion of fine-grained sediment (Dellwig et al., 2000).Increased macrozoobenthos abundances with higher mud content have also been reported by other studies (Armonies & Hellwig-Armonies, 1987;van der Wal, Lambert, et al., 2017).Armonies and Hellwig-Armonies (1987) found that abundance increased significantly towards the shore when sediments were muddy in contrast to other sediment types, where abundance peaks at mean high tide or was not increasing.However, it is important to note that this pattern is not universally valid.Beukema (1976) found the highest species richness at intermediate mud content and tidal level in the Dutch Wadden Sea.

Bosco et al. (2022) investigated how different diversity metrics
and biomass changed along a sediment texture gradient and found changes for all parameters with mud content and sediment grain size.
There are some methodological differences between the study of Bosco et al. (2022) and our study: Simpson diversity was calculated based on biomass and not abundance, GAMs were used to investigate the relationship, and in our study, the mud content had a broader range.Several studies found that macrozoobenthos species have specific preferences for certain sediment types (Anderson, 2008), resulting in a unimodal response of species richness to sediment grain size, which can be explained by an overlap of species that prefer muddy or sandy habitats at intermediate level grain sizes (Bosco et al., 2022).It is important to keep in mind that next to environmental filtering there are also other abiotic processes that can influence the macrozoobenthos community structure, for example, the spatial heterogeneity of abiotic factors (Kraft et al., 2015).

| Impact of local biotic filters (Chl a, TOC) and seasonality on macrozoobenthos: Abundance, biomass and biodiversity patterns (H4)
We expected that sampling stations with high but not extreme resource availability in terms of Chl a and TOC would support a higher abundance, biomass and diversity of macrozoobenthos, as microphytobenthos and organic carbon are the major food sources of macrozoobenthos (Christianen et al., 2017;Daggers et al., 2020;Riekenberg et al., 2022).We cannot confirm part of H4, as we could not observe any relationship of the community characteristics with Chl a (Table 2).Following the SAB model of Pearson and Rosenberg (1978), we expected that extreme values of TOC would lead to a decline in abundance, biomass, species richness and diversity, which we can confirm for all parameters except for abundance (Figure 3g-i, Table 2).Abundance is expected to peak shortly before extreme enrichment, as more opportunistic, small-sized species become abundant (Pearson & Rosenberg, 1978;Van Colen et al., 2012), and this might be the reason we did not observe any relationship between abundance and resource availability.
A study from Drent (2010) also could not link summer Chl a levels from the water column to biomass, abundance and diversity of macrozoobenthos.Unlike the Drent study which used phytoplankton as resource proxy, we used microphytobenthos for resource availability.Despite microphytobenthos being the most important food source for many macrozoobenthos species (except for filter feeders and predators; Christianen et al., 2017;Daggers et al., 2020;Lange et al., 2018), we found no significant effects.Furthermore, food resources are not always the limiting factor at all elevations.
At high elevations where the biomass of macrozoobenthos is low, the corresponding food demand is also low.Consequently, no food shortage is present, which makes other environmental factors more important (Beukema et al., 2002).It is also plausible that few superior competitive macrozoobenthos species dominate at the most productive environments, resulting in the observed negative relationship between macrozoobenthos species richness and productivity (Kondoh, 2001).
Additionally, factors such as the actual consumption rates of macrozoobenthos and the productivity of macrozoobenthos and microphytobenthos would be of interest, since top-down control by grazing from meiofauna and macrozoobenthos can influence the Chl a biomass at a site (Daggers et al., 2020).Information from these factors would reveal not only the existent resource availability, but also the total resource availability without consumption by other trophic levels.Top-down control from macrozoobenthos on microphytobenthos can change the spatial patterns of microphytobenthos in intertidal mud flats as shown in field experiments (Weerman et al., 2011).
The decrease in microphytobenthos coincided with an increase in the abundance of macrozoobenthos (Weerman et al., 2011).
We can partly confirm H4 as we found increasing abundances from spring to summer (Figure 3k, Table 2).Spawning and larvae settlement in spring results in higher abundances in summer when young specimen had time to grow to a size that was retained in sieves with a mesh size of 0.5 mm (Beukema & de Vlas, 1989).We could not detect any changes in biomass, probably because freshly settled individuals are still small, contributing only minor amounts to the total biomass (Beukema, 1974).Other factors, such as individual weight gains (Beukema, 1974), migration or the vulnerability of species sensitive to high summer temperatures (Beukema & Dekker, 2020), could act as balancing forces in influencing changes in biomass.We have to partly reject H4, as we did not predict that diversity decreased during summer.The drop in diversity can be a consequence of unusually high summer temperatures in 2019 (Meier et al., 2020;Zielinski et al., 2019), which might have caused migration or die-off of some species sensitive to high temperatures (Beukema et al., 2009).
Bearing in mind that macrozoobenthos species themselves have important ecosystem functions like bioturbation, destabilization or stabilization or reef formation that can influence the sediment characteristics (Eriksson et al., 2010;Meysman et al., 2006;Reise, 2002), physical, biochemical properties and hydrodynamics of their surroundings.Additionally, it is still unclear whether biotic resource filters, such as microphytobenthos Chl a biomass and TOC, which are consumed by macrozoobenthos, are good proxies for food availability (Hagerthey et al., 2002;Pratt et al., 2015;Weerman et al., 2011).
Thus, it is difficult to disentangle the cause-and-effect dynamics of environmental filters and resource filters on macrozoobenthos communities.

| CON CLUS ION
In conclusion, our research revealed new insights on how benthic macroinvertebrate abundance, biomass, species richness and diversity change with regional as well as local abiotic and biotic environmental filters in an important boundary zone between land and sea.
In view of climate change and rising sea levels, environmental abiotic and biotic conditions are expected to change in back-barrier tidal flats in terms of sediment composition, resource availability, elevation and extreme weather events.Interestingly, environmental factors shaped not only diversity patterns and dynamics of intertidal benthic macrozoobenthos on back-barrier tidal flats in the Wadden Sea but also community composition.As abiotic factors often covary, but with contrasting impacts on the macrozoobenthos community, it is challenging to entangle which abiotic factors are influencing the macrozoobenthos community in which way, especially when relying on correlative data.The regional and local relationships of biomass, abundance and diversity we observed are similar to those found in previous studies on macrozoobenthos in the Wadden Sea.
Despite its protection status, the Wadden Sea exists close to human pressures.Negative impacts do not need to have their source inside the protected area, which is for example the case for nutrient inputs from rivers (Hill et al., 2021).Management and protection of tidal flats worldwide have to be improved and monitored to protect the biodiversity and unique species assemblage that can be found in these coastal ecosystems.
Non-metric multidimensional scaling (NMDS) ordination diagram of the sampling stations (n = 162) at the three islands shown as colour coded (Nor = Norderney, Sp = Spiekeroog, Wang = Wangerooge).Effects of measured environmental variables on species assemblages of different East Frisian islands.Environmental drivers are indicated with vectors.The closer data points are to the centre the more similar they are in their species assemblage and the less they are influenced by the environmental parameters.Stress: 0.21, Dimension = 2.
), resembling the geographical distance of the islands.Some samples from Spiekeroog and Wangerooge clustered further apart from the rest of the samples.This part of the community was associated closely with the environmental factors mud, Chl a and TOC and the species Manayunkia aestuarina, Enchytraeidae, Pygospio elegans and Baltidrilus costatus.Other species such as Scoloplos armiger clustered on the opposite side of the mud vector and preferred sandy sites.
Results of models on biomass, species richness, abundance and inverse Simpson diversity over all islands and seasons in relation to environmental parameters (longitude, elevation, mud content, TOC, Chl a and season).
The marginal R 2 only considers the variance of the fixed effects, so it stands for the proportion of variance explained by the model's fixed effects, while the conditional R 2 lands have been, to varying extents, protected with dikes, groynes and wooden barriers.While on Spiekeroog no coastal protection has been established, Norderney is protected with wooden barriers and Wangerooge is protected with groynes (NLWKN, 2010).The separation of salt marsh and mud flat by different flood protection measures might have contributed to the differences in community composition among the islands.Furthermore, local abiotic factors can also change with longitude, like the grain size and the mud content, and can have an influence on the community composition.Compton et al. (2013) found changes in the macrozoobenthos assemblages from west to east in the Dutch Wadden Results of PERMANOVA (adonis2), for macrozoobenthos relative abundances (Bray-Curtis dissimilarities, 9999 permutations) to investigate the effect of island on the macrozoobenthos community composition on the intertidal flats of the German Wadden Sea.
TA B L E 1 Mean abundance and standard deviation and mean biomass and standard deviation per island and species per core (0.0078 m 2 ).